Zigzag graphene nanoribbon edge reconstruction with Stone- Wales defects 
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In this article, we study zigzag graphene nanoribbons with edges reconstructed with Stone- Wales 
defects, by means of an empirical (first-neighbor) tight-binding method, with parameters determined 
by ab-initio calculations of very narrow ribbons. We explore the characteristics of the electronic 
band structure with a focus on the nature of edge states. Edge reconstruction allows the appearance 
of a new type of egde states. They are dispersive, with non-zero amplitudes in both sub-lattices; 
furthermore, the amplitudes have two components that decrease with different decay lengths with 
the distance from the edge; at the Dirac points one of these lengths diverges, whereas the other 
remains finite, of the order of the lattice parameter. We trace this curious effect to the doubling of 
the unit cell along the edge, brought about by the edge reconstruction. In the presence of a magnetic 
field, the zero-energy Landau level is no longer degenerate with edge states as in the case of pristine 
zigzag ribbon. 

PACS numbers: 81.05.ue,72.80.Vp,78.67.Wj 
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I. INTRODUCTION 

At the present time, the most promising scalable 
growth methods of graphene films are either based on 
epitaxial growth on silicon carbide^ ^ or on chemical va- 
por deposition (CVD) of graphene on metal surfaces.!^^ 
Yet, the latter methods do not produce graphene films 
with electronic mobilities as high as those reported for 
exfoliated graphene.^ ^ Electronic transport^ ^^ in CVD- 
grown graphene is hindered by grains, grain boundaries 
and atomic patchwork quilts, ^^ ^^ which can be inter- 
preted as topological defects A^^ 

CVD-grown materials are in general polycrystalline in 
nature, having their physical properties dominated by 
the grain boundaries' size. The situation is no differ- 
ent for graphene. For this material, it is theoretically 
expected that some of its electronic properties will be 
markedly different from its exfoliated counterpart, as sug- 
gested by calculations of formation energies of different 
types of grain boundarie^^ and by the transport mea- 
surements and theoretical calculation^^ in high-quality 
CVD-growiP graphene. 

Due to graphene's hexagonal structure, originated 
from the sp^ bonds, the grain boundaries are expected to 
be formed of pentagons- heptagons pairs, known as Stone- 
Wales (SW) defects. ^'^ Recent atomic resolution TEM 
studies^^ ^^ ^^ have allowed to visualize grain boundaries 
in CVD-grown graphene. These experimental studies 
have shown that the grain boundaries are not perfectly 
straight lines and that the 5-7 defects along the bound- 
aries are not periodic. These type of defects have a pro- 
found effect on the threshold for mechanical failure of 
the graphene membranes, which is reduced by an order 
of magnitude, relatively to the exfoliated membranes. 
In what concerns the electronic properties, it has been 



shown that the measured electronic m obilities depend on 
the details of the CVD-growth recipes .I^EI^El 

Furthermore, as shown by recent TEM studies,"^ these 
extended pentagons-heptagons pairs defect lines inter- 
cept each other at random angles, forming irregular poly- 
gons with edges showing stochastic distribution of length, 
making it extremely difficult to make theoretical studies 
of these defects using microscopic tight-binding models. 
On a different tone, it has been argued that these de- 
fect l ines c an act as one-dimensional conducting charged 
wires.^^^^ The charging of these topological wires is 
achieved by the self-doping mechanism. ^^ 

As said, studying Stone- Wales defects in the bulk of 
graphene, using microscopic tight-binding models, is a 
difficult task, due to the breaking of translational geom- 
etry. On the other hand, the grain boundary formed 
by the 5-7 defect lines effectively cre ate an edge, giving 
rise to an enhanced density of state^*^ at the Dirac 
point, all in all equal to what is found at the edges of 
zigzag nanoribbons.^^ ^^ Evading the difficulty of study- 
ing topological defects in the bulk of graphene, we take, 
in this article, the approach of studying the formation of 
Stone- Wales defects at the edges of zigzag nanoribbons, 
supported by the experimental findings that grain bound- 
aries effectively act as edges of the crystalline grain.l^^^ 
We will be focusing our study on the electronic properties 
of graphene nanoribbons close to the Dirac point, for the 
effect of Stone- Wales defects have their largest impact on 
the properties of graphene at low energies. 

Ab-initio calculations have shown that when SW 
defects are present in graphene nanoribbons (GNRs), 
the energy decreases as the defect gets closer to the 
edge of the ribbon.^^ Other first principles studies have 
shown that the formation of SW defects at the edges 
of both armchair and zigzag GNRs (respectively, AG- 



NRs and ZGNRs), stabilize them both energeticahy 
and mechanicanyP^Ji^ The zigzag edge, in particular, is 
metastable under total reconstruction with SW defects, 
and a planar reconstruction spontaneously takes place at 
room temperature.^^ ^^ 

Edge-reconstructed ZGNRs by means of SW defects, 
are claimed to be stable only at very low hydrogen pres- 
sure (well below the hydrogen pressure at ambient con- 
ditions) and very low temperatures!^ However, recon- 
structions of the zigzag (as well as armchair) edges have 
been recently observed with high-resolution TEMP^H^ 
The recent work of Suenaga et al.^ on single-atom 
spectroscopy using low- voltage STEM, may be used as 
yet another mean of identifying edge reconstructions of 
graphene ribbons. Moreover, refinements in other tech- 
niques, such as Raman spectra of the edges, ^^ STM im- 
ages of the edges,'^and coherent electron focusing,!^ may 
help in the identification of these kinds of edge recon- 
structions. 

In this work, we have studied various reconstructions 
of zigzag edges with SW defects, namely zz(57), zz(576) 
and zz(5766) (see Fig. fl]). However, in this article, we 
give special emphasis to the case of total reconstruction of 
the zigzag edges, zz{b7), because it is the most sta ble con- 
figuration in the absence of hydrogen passivation! ^^ * ^^ * ^ -^ 
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Figure 1: SW-reconstructed zigzag edges: the pristine zigzag 
edge in (a), also known as zz edge; the zz(57) edge in (b); the 
zz{576) in (c); the zz{5766) in (d); 

This article is organized as follows: In section [Hj we 
study the electronic structure of wide zigzag ribbons, 
with edges reconstructed by Stone- Wales defects, using 
an empirical tight-binding model. In Subsection |II A[ 
we start by computing the model parameters using the 
results of ab-initio simulations. Based on the empiri- 



cal tight-binding model presented in Subsection |IIB[ we 
study the electronic structure of ZGNRs whose edges 
were reconstructed by SW defects, with a focus on the 
edge states showing up in these systems (Subsection 
II C ) . We find that some modifications (relatively to the 
pristine ZGNR) are introduced in the electronic struc- 
ture, as well as, that the edge states of the edge recon- 
structed ZGNRs are distinct from those of the pristine 
ZGNRs. Finally, in Subsection II D[ we explore the im- 



plications of the presence of a magnetic field directed per- 
pendicularly to the ribbon plane, in the electronic struc- 
ture and edge states of a zigzag ribbon, pinpointing the 
modifications originating from the edge reconstruction. 



II. TIGHT-BINDING STUDY OF RIBBONS 
WITH RECONSTRUCTED EDGES 

In this Section, the main issue under discussion is 
the behavior of the edge states of wide zigzag ribbons, 
whose edges have been reconstructed due to the forma- 
tion of edge Stone- Wales defects (see Fig. [l]). The huge 
amount of computational resources needed to study large 
physical systems employing ab-initio methods, makes it 
prohibitive to explore the physics of wide ribbons using 
such techniques. An alternative approach, is to use phe- 
nomenological tight-binding models, which being com- 
putationally not so demanding also give a microscopic 
understanding of the electronic properties of these types 
of systems. 

In order to provide an accurate (tight-binding) descrip- 
tion of the reconstructed edges, we start by performing 
ab-initio simulations of narrow ribbons, from which we 
extract the values of the hopping amplitudes at the edge. 
These hopping amplitudes are posteriorly used in the 
construction of the tight-binding Hamiltonian, in which 
the study of the edge states for large ribbons is based on. 



A. Parametrization of the hopping amplitudes 
using ab-initio methods 

We used Density Functional Theory (DFT) to 
parametrize the tight-binding. The calculations were per- 
formed using the code AIMPRO,^ under the Local Den- 
sity Approximation. The Brillouin-zone (BZ) was sam- 
pled for integrations according to the scheme proposed by 
Monkhorst-Pack.^^ The core states were accounted for 
by using the dual-space separable pseudo-potentials by 
Hartwigsen, Goedecker, and Hutter.'^The valence states 
were expanded over a set of s-, p-, and d-like Cartesian- 
Gaussian Bloch atom-centered functions. The k-point 
sampling was 16 x 2 x 1 and the atoms were relaxed 
in order to find the equilibrium positions. A supercell 
with orthorhombic symmetry was used; the cell parame- 
ter in the infinite direction was 4.885^4. A vacuum layer 
of 12.7 A in the nanoribbon plane and 10. GA in the nor- 
mal direction were used in order to avoid interactions 
between nanoribbons in different unit cells. 

In what follows, we will focus on the ZGNR with to- 
tally reconstructed edges, the most stable of this family 
of reconstr uctions i n the absence of hydrogen passivation 
(see Fig. [iJ)J 26 | 27 | 3Q ]^q^^ ^^^^ ^^^ dangling bonds that are 
on the origin of the zigzag edge react iveness, are elimi- 
nated by the reconstruction of the edge, forming triple 
bonds between the outer carbon atoms at the edges (/12 
bond in Fig. [4]). In the literature, the SW totally recon- 



structed edge is usually named as zz{b7). Note that the 
unit cell of such a ZGNR has twice the size of the unit 
cell of the pristine ZGNR (see Fig. [3|. The generaliza- 
tion of the following study for SW edge reconstructions, 
other than zz{57), for example zz(576), zz(5766), etc., is 
straightforward. 

In Fig. [2J we show the relaxed edge geometry of a to- 
tally reconstructed edge (in absence of hydrogen passiva- 
tion), zz{57), as obtained from the ab-initio calculations, 
together with the inter-carbon distances and the angles 
between carbon bonds. 




Figure 2: Relaxed edge geometry of the totally reconstructed 
zigzag edge (in absence of hydrogen passivation), zz{57), ob- 
tained from the ab-initio calculations. The numbers refer to 
the bond lengths in angstroms; the capital letters refer to the 
angles between two adjacent bonds. 

The procedure for determining the hopping amplitudes 
at the reconstructed edge is the following. From ab-initio 
calculations one obtains the different carbon-carbon dis- 
tances at the edges of the ribbons as well as the values of 
the angles between carbon bonds (see Fig. [2]) . In the case 
of the zz{57) edge reconstruction, our first principles cal- 
culations show that the ribbons remain planar (we have 
allowed the system to relax along the three spatial di- 
mensions); therefore the values of the angles in Fig. |2] 
play no role in the determination of the hopping ampli- 
tudes, since these arise from ppn hybridization. Using 
the carbon-carbon distances, we compute the hopping 
amplitudes using the parametrization^ 
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where r^j stands for the distance between the carbons la- 
beled by i and j (given in unities of Angstroms), the 
adimensional parameters 0^2 = 1.2785, 0^3 = 0.1383, 
Qf4 = 3.4490, while a^ is the carbon-carbon distance in 
the bulk (in units of Angstroms). ^^ The hopping renor- 
malizations, /i^, h[ and Vi^ (see Fig. |4]for their definition) 
are given by the r{rij) for the corresponding carbon- 
carbon distances at the edge. For the zz{57) edge, in 
the absence of passivation, the values of these renormal- 
izations are listed in Table [H 



B. The Tight-Binding Hamiltonian of a ZGNR 
with zz{57) edges 

The simplest model one can construct describing non- 
interacting electrons in a ZGNR whose edges have been 
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Table I: Values of the hoppings in unities of t (which we also 
call hopping renormalizations) for a zz{b7) (see Fig. Inland 
Fig. B calculated from the C-C distances obtained from the 
DFT numerical calculations using Eq. (IT]). 



reconstructed by SW defects is the first neighbor tight- 
binding (TB) model. Generically, a ribbon has TV zigzag 
rows of atoms along the longitudinal direction (0 < n < 
N — 1) and, in each unit cell, there are P zigzag columns 
of atoms {I < p < P). In the case of a zz{57) edge, 
P = 2. 

The TB Hamiltonian for the edge reconstructed 
ZGNR, can be written as 
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where H^ stands for the Hamiltonian of the region in the 
vicinity of the upper edge of the ribbon (at n = in Fig. 
p|, H^ stands for the region in the vicinity of the lower 
edge (at n = A/" — 1 in Fig. [3| and H^'^^^ stands for the 
bulk of the ribbon. 
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Figure 3: Scheme of a ZGNR with its edges totally recon- 
structed by SW defects [a zz{57) ribbon]. Details of the edge 
are represented in Fig. [4] 

The ab-initio results (see Fig. |2]) show that only in the 
first two rows are the hopping parameters between two 
adjacent carbon atoms different from their usual value 
in the bulk. Thus, we choose to identify term H^ {H^) 
in the full Hamiltonian with the two upper (lower) rows 
of atoms of the ribbon. The annihilation operators of 
the four numbered atoms in row n = (see Fig. [4|, are 
denoted by di{m)^ d2{m)^ ds{m)^ and ^4(771), while those 
referring to the four numbered atoms in row n = 1, are 
denoted by ci(m), C2(m), cs{m), and 04(171). 

For the sake of clarity, we will separate in H^ the 
terms referring to each row, n = and n = 1, and to the 
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Figure 4: Details of the ZGNR with its edges totally recon- 
structed by SW defects, zz{57). The hi, hi and Vi stand for 
the factors giving the renormalization of the hoppings, t, be- 
tween nearest neighbors in the vicinity of the defect. We omit 
the lower edge because it is analogous to the upper one. 



coupling between them, H^ = H^^q + H^^^ + H^ . For 
row n = 0, we have 



h: 



n=0 



m ^ i=l 



{m)di-^i{m) 



-\- h4d[{rn)di{rn -\- 1) -\-h.c.>, 



(3) 



while for row n = 1, 



H^^, = -i^|^[ft^cl(m)c,+i(m)' 

m ^ i=l 

+ h'^cl{m)ci{m ^ 1) + h.c. L 



(4) 



and for the coupling between row n = and row n = 1, 



H^ = —ty ^ \vic\{m,)di{m) + V2c\{m)di{m) 



h.c. 

(5) 



Recall from Table |l] that hi = /13, h[ = h'2^ /13 = /14, and 

Vi = V2- 

The term H^'^^^^ corresponding to the Hamiltonian of 



the bulk (between row n = 2 and n = A^ — 3), is given by 
Ar-3 



i7' 



bulk 



m n—2 ^ 

+ aj(r?i; n + 1) 61(771; n) + a2(m;n) 

+ a|(7Ti + 1; n) + a2(m; n + 1) h2{m]n) ^\i.c.\^ 



(6) 

where ap{m;n) [bp{m;n)] is the annihilation operator of 
an electron state localized at the atom of sub-lattice A 
(B) in column p (p = 1, 2 for a zz{57) edge) and row n, 
of the unit cell labeled by m. The term H^ describing 
the lower edge is analogous to the upper edge term, H^ . 
Recall that the /i^, v and h[ parameters in the equations 
defining the tight-binding Hamiltonian correspond to the 
values of the hoppings in units of t. In addition, we make 
the following identifications: 

(ii(4)(m) -^ 6i(2)(^;0), 
d2i3){m) -^ ai(2)(^;0), 
ci(3)(^) -^ «l(2)(^;l), 

C2(4)(^) -^ ^1(2) (^;1)- 

With no loss of generality, we assume periodic bound- 
ary conditions along the ribbon x-direction. This sim- 
plification allows us to diagonalize the Hamiltonian with 
respect to m by Fourier transforming H along the x- 
direction, 



^ = E^'= = E[^^ 
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bulk 
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Having determined the values of the hi , v and h[ (see 
Table |l]) , we compare in Fig. [s] the obtained low-energy 
spectrum from the ab-initio calculations with that re- 
sulting from the numerical diagonalization^^ of the tight- 
binding Hamiltonian H^^ Eq. ([t]). 

As we can see in Fig. [5J the DFT and TB numer- 
ical calculations for narrow zigzag ribbons with zz{57) 
edges originate low-energy spectra with similar features. 
The differences between the DFT and the TB spectra 
are probably due both to the simplified character of the 
TB treatment (especially the first- neighbor approxima- 
tion) and to finite size effects affecting both systems dif- 
ferently. In fact, it is well known that even for an ac- 
curate description of ab-initio of bulk graphene bands, 
one needs a tight-binding model including hoppings up to 
third-nearest neighbors. ^^ Since our interest is the under- 
standing of the main features of the low-energy spectra, 
we keep in the tight-binding model only the first neighbor 
hopping. 

In the reduced Brillouin zone, arising from the dou- 
bling of the unit cell along the edge (x) direction, 
the Dirac points of bulk graphene appear at K = 
27r(l/2, - V3/2)/3 and K' = 27r(-l/2, V3/2)/3. In a rib- 
bon, they will show up at ka = 7r/3 and at ka = — 7r/3, 
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Figure 5: Comparison between the low-energy spectrum of a 
ZGNR (with both edges reconstructed) obtained with DFT 
(blue) and TB (red). In (a) the ribbon has a width of 18A 
(or N — % zigzag lines), while in (6) the ribbon has a width of 
31A (or A^ = 14 zigzag lines). The Fermi level is at Ejt = 0. 



where k is the momentum along the edge direction. 
We now focus on the dispersive energy levels present 
around the Fermi level, appearing between these two 
Dirac points. 



C. Edge states of a zzi^l^ edge 

In a finite graphene sheet, energy levels appearing 
outside the range of allowed electronic states of bulk 
graphene correspond to states localized at the edges, 
called 'edge states', as usual in surface physics. Conse- 
quently, from Fig. |6|a), we can guess that zz{57) edges 
allow both high and low-energy edge states (respectively, 
the levels h and / in Fig. [6]); this contrasts with what 
happens in the pristine zigzag edge (only low-energy edge 
states) ) ^^tt24| jj^ what follows, we will focus on the physi- 
cally relevant low-energy ones. 




Figure 6: (a): Low-energy spectrum of a ZGNR with N — 
30 and both edges totally reconstructed with SW defects, 
zz(57). Remember that the energy spectrum of a totally re- 
constructed edge, results from a doubled unit cell relatively to 
the pristine ZGNRs, and consequently is double- folded rela- 
tively to the latter. The labeled energy levels (b): Low-energy 
spectrum of a pristine ZGNR with A^ = 30 and zigzag edges. 
Both (a) and (b) were obtained using a tight-binding model. 
In both cases, the shaded areas indicate which levels are al- 
lowed in bulk graphene. Eigenstates corresponding to levels 
that are outside the shaded region, will necessarily be located 
at the edges of the ribbon both in (a) zz{57) edges and in (b) 
perfect zigzag edges. In the (a) panel, the levels labeled by h 
and /, stand, respectively, for high and low-energy. 



Since edge states decay exponentially into the bulk, 
in wide ribbons they can be studied as states of semi- 
infinite ribbons: states at different edges are uncoupled 
if the ribbon width is much larger than the decay length. 



In the case of a semi-infinite ribbon with pristi ne zigzag 
edges, the edge states occur at zero-energy !^^*^^* In such a 
case, the tight-binding equations simplify to independent 
recurrence relations for the amplitudes of the A and B 
sub-lattices, which yield, transparently, the exact wave- 
functions, the analytical expression of the decay length as 
function of the momentum along the edge, and the range 
of momentum values in which such states are possible. 
In the case of a zigzag ribbon with edges totally recon- 
structed with SW defects, we face the complication that 
edge states have dispersion^ and are not at zero energy. 

To investigate the possibility of low energy edge states 
of such a system, we must solve the Schrodinger equation. 



Hk\li,k) 



^^,k 



\li,k) 



(8) 



for \e^^k\ /^ <^ 1, where Hk is the same as that obtained 
from the transformation in Eq. [7J of the Hamiltonian 
given by Eqs. ([2])-([6| with n > 0. Note that Hk defines 
effectively a ID problem in the transverse direction of 
the ribbon. Consequently, we can express any eigenstate 
l/i, k) as a linear combination of the site amplitudes along 
the transverse direction of the ribbon, 



N-l 2 



1^'^) = X^ X] [^p(^;^)l«;^; 



p,n) 



n=0 p=l 

Bp{k;n)\b;k;p,n) , 



(9) 



with the one-particle states, |r;/c;p,n) = r^(/c;n)|0), 
where p = 1, 2 and n = 1, . . . , A/" — 1 define, respectively. 



the column and the hne of the unit ceh, and r = a^b. To 
Hghten our notation, we have identified the states at the 
upper edge as 



Mi(4);fc) 


= |5;fc;l(2),l), 


(10a) 


M2(3);fc) 


= |a;fc;l(2),0), 


(10b) 


ki(3);fc) 


= |a;fc;l(2),l>, 


(10c) 


|c2(4);fc) 


= |6;fc;l(2),2), 


(lOd) 


while at the lower edge 






|C5(7);fc) = 


\a;k;l{2),N-2), 


(11a) 


|C6(8);fc) = 


\b;k;l{2),N-2), 


(lib) 


M5(8);fc) = 


\a;k;l{2),N-l), 


(lie) 


14(7); fc) = 


\b;k;l{2),N-l). 


(lid) 



Note that there are four states per zigzag row (identified 
by n), coming from the four sub-lattices Ai, 5i, A2 and 
B2. Equating coefficients, we obtain a set of 2 x 2 x 
N (tight-binding) equations, where N is the number of 
zigzag rows of atoms in the unit cell. 

To build an analytical description for edge states 
in a semi-infinite ribbon, with row index n > 0, 
we write the TB equations in matrix form, where 
A{k;n) = [Ai{k]n),A2{k;n)] and B{k;n) = 

[Bi{k; n), B2{k] n)] will stand for column vectors. 

For rows with n > 1, the relations between the ampli- 
tudes are the same as those of a pristine ribbon: 



A(A;;n + l) 
B{k;n) 



WAA{k;n) 
WBB{k;n- 



= -(^)B(^;n + l), (12a) 
l) = -(^)A(fc;n). (12b) 



The matrices Wa and Wb, defined in Eqs. ( |A11| ), 
commute and, therefore, share a common eigenbasis, 
{n+,n~} (see Appendix [A| for details). Let us denote 
the corresponding eigenvalues by ^^ and ^^ , respectively. 
These quantities depend on the value of the longitudinal 
momentum, k and are given by: 






ts = 

G = 



-2cos{ka/2)e'^, (13a) 

2ism{ka/2)e'^, (13b) 

-2cos(te/2)e-^^ = (^+)* , (13c) 

-2isin(te/2)e-^* = (^^)* . (13d) 

Changing to the {u~^^u~} basis, 

A{k;n) = a-^{k;n)u~^ -\- a-{k;n)u~ ^ (14a) 

B{k;n) = /3+(A:; n)n+ + /3_(A:; n)^-, (14b) 

we obtain Eqs. (IT2| in the form. 



aa{k; n + 1) - G<^a(^; n) = - (^- j (3^{k; n + 1), 



(15a) 



P.{k', n) - {arm n + 1) = - (^) a^{k; n), (15b) 



where a = ±1. Note that with the two possible values 
for cr, Eqs. (15) give four equations. These equations 



describe two independent ID AB chains in the n coordi- 
nate, one for each of the modes u+ and u_; the hopping 
amplitude alternates between —t and t^^. 

The two modes u+ and u_ are easily interpreted. If 
we look for propagating solutions {qct real). 



aa{k;n) = a^(/c)e'^"'^, 
/3,(fc;n)=/3,(^)e^^-^ 



we arrive at the equation 
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(16a) 
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(17) 



Low energy states correspond to (e/t)^ <C 1; but it can 
easily be checked from Eqs. (13), that |^^| > ^/2^ for 



all ka in the F.B.Z., whereas^^ | ^ 1 around ka = 
±7r/3. Hence, propagating states of the a = -\- modes 
have an energy of order t; the a = — modes are the low 
energy bulk states when k is near the Dirac value. The 
existence of these two modes reffects the folding of the 
Brillouin zone to account for the doubling of the unit cell. 
At the Bloch momentum of a Dirac point there are two 
different energy levels, only one of which is of low energy, 
and corresponds to the u_ mode. In fact, inspecting 
the relation between the Ai and A2 a mplitudes in the 
u_mode [see Appendix \A\ Eqs. (A17)] one sees that it 
corresponds to what is expected from a plane wave at a 
Dirac point. 

Nevertheless, for decaying states {qct with an imaginary 
part), we cannot exclude the possibility that low energy 
states can have a a = + component, because in that 



case, the right hand side of Eq. (17) has a factor (1 



-^q+ 



(e 



aT 



mq+ 



), which can be close to zero. We will see 



in a moment that the boundary conditions (BCs) arising 

from the zz{57) edge bring about precisely this situation. 

Let us now discuss what kind of solutions are obtained 



from Eqs. (15) if the system supports zero energy states. 



For zero energy, the bulk Eqs. (15) become independent 
recursion relations 



acr(/c;n + l) = aa(j{k;n 
1 



/3^(/c;n + l) 



i^A 



^/3a{k;n). 



(18a) 
(18b) 



From Eqs. (13), l^^l > >/2, thus requiring Qf+(A:,n) = 0, 



otherwise we would have a non-normalizable state. Also, 
we must have either a_(/c, n) or /3_(/c, n) = 0, depending 
on whether |^^ | is greater or smaller than 1. Consider, for 
instance, the latter case: the required conditions for zero 
energy states would then be a+(/c, n) = /3_(/c, n) = 0. 

The previous paragraph did not impose any type of 
conditions arising from the boundary. It turns out that 
the existence of zero energy states depends on the specific 
form of the boundary conditions. We note however, that 
in this type of edge reconstruction surface states always 
exist, but not necessarily at zero energy. The boundary 



conditions can be derived form the tight-binding equa- 
tions for the rows n = 0, 1. As shown in Appendix [A} 
Eq. (A19), they can be approximated by the zero energy 
BCs, cx{k; 2) = A^/3(/c;2), where M is Si k dependent 
matrix defined exphcitly in Appendix [Aj in fuh, 

a^{k; 2) = >(++/3+(/c; 2) + M+_/3_(/c; 2), (19a) 
a-{k; 2) = M-^(3^{k; 2) + M—(3-{k; 2). (19b) 

In the case where zero energy states exist, the boundary 
conditions defined by Eqs. (fT9l) are exact. For the k 



values for which |^^| < 1, zero energy states require, 
as we have seen, a+(A:) = f3-{k) = 0; this is possible 
only if A^++ = 0. This condition is, in fact, verified in 
certain limits, the simplest one corresponding to ignoring 
the hopping renormalizations at the edge, that is, taking 
hi = V = h^, = I'm which case the matrix M reads 



M = -4sin^(/ca) 



(G) 



Another interesting limit to consider is h[ 
case, one obtains 



(20) 



1. In this 



(21) 



and consequently, zero-energy states should be observed 
if hi - h2h^ = 0. 

We have confirmed these results by numerical diag- 
onalization of the tight-binding HamiltonianP^ In both 
situations, as A^++ = 0, the zero-energy states appear 
in the range where |^^| < 1, i.e., \ka\ < 7r/3, and have 
the form (for n > 1) 

a-{k]n) = 

pj,{k',n) = 




with 



a_{k) = -4sin2(te)(^+)*/3_(A:). 



(22a) 
(22b) 

(23) 



In Fig. [7|we compare numerical diagonalization results 
with those of the present analysis, for the simplified sit- 
uation where hopping renormalizations at the edge are 
ignored, hi = v = h[ = 1.^^ The squared amplitudes of 
the edge states, of a narrow ribbon with A^ = 30 (65A 
wide), calculated numerically, are indeed in very good 
agreement with those of the edge states of a semi-infinite 
ribbon obtained analytically, from Eqs. (22) and (23). 



Unlike the zero energy states occurring in unrecon- 
structed ZGNR, the wave function amplitudes of the edge 
states are non-zero in both sub-lattices. Those familiar 
with the Dirac equation description of graphene might 
find this result surprising, since, at zero energy, the equa- 
tions for the A and B fields decouple, and only one of 
them can be non-zero. ^^ However, as can be seen in Fig. 
u\ panels (a4)-(a5) - which refer to a value of k close to a 



Dirac point -, the decay length is much shorter in the B 
sub-lattice; this is related to the fact that the B ampli- 
tudes correspond to the a = + mode, which, in the bulk, 
is high energy, and has a finite decay length, of the order 
of a single row width, even at the Dirac point, contrast- 
ing with the (J = — mode, whose decay length diverges at 
the Dirac point. So, away from the boundary, the edge 
state wave function is, in fact, similar to that of a ZGNR, 
because the amplitude at the B sub-lattice is exponen- 
tially smaller than in the A one; but the reconstructed 
edge requires the presence of the confined a = + mode, 
in order to satisfy the BC. When we move away from the 
Dirac point, the distinction between high and low energy 
modes washes away, and both modes are confined within 
atomic distances to the edges [Fig. [tI panels (a2)-(a3)]. 

At this point we come back to the consideration of real 
edges, where the hopping parameters have the values in 
Table fl| One does not find A1++ = 0, and the BCs of 
Eqs. (119]) are no longer compatible with the conditions 
for zero energy states, Qf+(/c) = P-{k) = 0, (or a+(/c) = 
a-{k) = 0, if l^^l > 1); edge states, if they exist, have 
to be dispersive. The dispersiveness of the edge states' 
levels can be seen in panel (61) of Fig. ItI 

We analyze this situation in detail in Appendix [BJ If a 
semi-infinite ID AB chain has a zero energy edge state 
with BC, say B{0) =0, it will still have a low but finite 
energy one, if the BC is replaced by B{0) = sA{0) with 
1^1 < 1. In the present case the situation is more complex, 
because the problem involves two ID chains [Eqs. ([iS])], 
coupled by the BC [Eq. (p!9|)]. The main conclusion still 
holds, and we expect low energy, dispersive, edge states 
near the Dirac points {ka = ±7r/3). 

In Fig. [S] we compare analytical results for a semi- 
infinite chain, obtained with the procedure described in 
Appendix [BJ with numerical diagonalization of a very 
wide ribbon {N = 600). The use of the zero energy 
BC of Eq. ([19]) correctly accounts for the wave func- 
tion and for the energy dispersion as a function of /c, 
but only very close to the Dirac point. It quickly devi- 
ates strongly from the numerical results as we move away 
from the Dirac point. This is to be expected, not only 
as a result of the violation of the low energy condition, 
but, more importantly, because the description in terms 
of bulk equations and simplified BCs will not hold when 
the edge state has such a short decay length that it lives 
mostly at the edge. Moreover, as stated before, near the 
Dirac points the localization length of the mode a = — 
diverges. As a consequence, the analytical analysis de- 
veloped in Appendix [BJ will only accurately describe the 
physics of zz{57) edged ribbons near the Dirac points if 
the ribbons are large. 

We can summarize the results of this subsection, saying 
that, as a consequence of the duplication of the unit cell. 
Stone- Wales reconstructed edges present a new type of 
edge state. 



A{k;n) = e^^+^^-'^a+u+ + a_e^^-^^-'^u_, (24) 



B{k;n) = e^^+^^-2)/3+u+ +/3+e^^-(^-2)u_, (25) 




0.15 



0.02 



Figure 7: (a): ZGNR with two simplified zz{57) edges {hi = 
V = hi = 1) and a width of 65 A (or A^ = 30 zigzag lines). The 
panel (al) shows the tight-binding low-energy spectrum in 
the FBZ. The two lowest-energy levels are colored in blue and 
red. The dashed (orange) horizontal line, signals the position 
of the Fermi level. Panels (a2) and (a4), show the edge state 
squared amplitude corresponding to the blue level in (al), 
for ka — 0.704 and ka — 0.842 respectively (whose position is 
identified in panel (a) by the vertical dashed green lines). The 
continuous (dashed) dark blue line stands for the amplitude in 
the Al -sub-lattice (Bi -sub-lattice) corresponding to the blue 
level in (al) obtained from the numerical diagonalization of 
the TB Hamiltonian (the red level is an identical edge state). 
Only the amplitudes Ai{k]n) and Bi{k;n) were plotted, be- 
cause A2{k;n) and B2{k;n) are identical to the former. The 
continuous (dashed) light blue line stands for the zero-energy 
edge state amplitude in the Ai -sub-lattice (Bi -sub-lattice) 
obtained analytically in a semi-infinite ribbon. Note the ex- 
treme coincidence between the numerical and the analytical 
edge states. Panels (a3) and (a5), show the same plots as in 
(a2) and (a4), but now with logarithmic scale in the ^-axis, to 
display the exponential decay of the squared amplitudes, (b): 
ZGNR with two real zz{57) edges (see Table IT]) and a width 
of 6b A. The panels are organized as those of (a). 




Figure 8: Comparison between the edge state levels obtained 
from numerical diagonalizing the tight-binding Hamiltonian 
of a ribbon with N = 600 {^ 900 A wide) (in blue), and the 
edge state level resulting from analytically solving the TB 
equations, by using simplified (zero-energy) BCs (in red), (a) 
shows the energy as function of k; panels (bl) and (c), show 
the square of the amplitudes as function of distance to the 
edge, for ka = —1.005 and ka = —0.804 respectively. These 
values of ka are identified in panel (a) by vertical green dashed 
lines; (62) panel is the same as (bl), but now in logarithmic 
scale. 



with the following features: (i) the states are dispersive; 
(ii) the wave- function, even for the semi-infinite ribbon, 
has non-zero amplitudes on both sub-lattices; (iii) close 
to k = ±7r/3, the Dirac points, the wave function ampli- 
tudes have two components decaying with very different 
rates, '^q- and ^g+, the latter remaining finite even at 
the Dirac point, and corresponding to a mode with only 
atomic scale penetration into the bulk. 

This last feature is strikingly apparent in Fig. [8J panel 
(62), where the faster decaying component in the B lat- 
tice dominates the wave function close to the edge, be- 
cause of a larger initial amplitude, |/3+| :^ |/3_|, but is 
supplanted by the one with slower decay, around n ~ 10. 



D. Perpendicular magnetic field 

When a perpendicular magnetic field is applied to a 
graphene sheet, electrons acquire a cyclotron motion, 
with quantized cyclotron radius and macroscopically de- 
generate energy levels, the so called Landau levels (LL). 
In a ribbon, LL degeneracy is partially lifted, because 
the edges interrupt the cyclotron orbits located close to 
them. In this section we discuss the effect of a perpen- 
dicular magnetic field in the low energy spectrum of the 
tight-binding models we have been discussing (a ribbon 
with a zz{b7) reconstruction). 

The introduction of a static magnetic field, applied per- 
pendicularly to the ribbon, B = Bcz, can be achieved 




and hexagons near the edge. Therefore, by Eq. (27), the 
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Figure 9: Peierls phases of a zigzag ribbon with totally recon- 
structed edges {N = 10). 



by a Peierls substitution of the hopping integrals,'^^'^ 

Uj ^ t,,e^2-^-, (26) 

where tij stands for the hopping integral between the 
position Ri and the position Rj in the absence of a mag- 
netic field, and the phase (j)ij is given by the line integral 



1 



XX o 



A ■ dr. 



^0 JR, 



(27) 



where A is the potential vector and (J)q = h/e is the flux 
quantum. Note that the magnetic flux through the area 
S, in units of the flux quantum ^o, is 



^ f da-B 

<Po is 



^(fdr-A= J2 ^'i- (28) 



Peierls phases around the edges are distinct from those in 
the ribbon bulk. We choose a gauge that yields Peierls' 
phases as shown in Fig. |9| clearly satisfying Eq. (28), (pQ 



being the magnetic fiux per hexagon in the bulk graphene 
lattice. 

The spectrum shown in Fig. floFa) is essentially the 
same as for a pristine ZGNR (apart from the folding of 
the Brillouin zone), the most prominent feature being a 
doubly degenerate zero energy level occurring between 
the two Dirac points. But what is displayed is, in fact, 
the spectrum of a ribbon with simplified zz{b7) edges, 
where hopping renormalizations were ignored {hi = 1, 
Vi = 1), and the pentagons and heptagons considered to 
have the same area as all the hexagons. 
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Figure 10: Panel (a) shows the tight-binding energy spectrum 
of a zigzag graphene nanoribbon with simplified 2:2; (57) edges 
(hi = V = hi = 1 and 0^ = 0^ = 20^6 = 206) with a 
width of 214A (or N = 100 zigzag rows) in the presence of 
a perpendicular magnetic field, B = SOT. The green dashed 
vertical lines in (a), indicate the different values of ka for 
which the edge states were plotted in {b)-{p). Panels {b)-{p) 
show in dark blue, for different values of ka, the wave function 
squared modulus of the two lowest-energy levels highlighted 
in panel (a) with blue and red fill. The light blue curves in 
panels (/)-(o) stand for the edge states obtained analytically 
for values of ka for which their energy is zero [see panel (a)]. 



The zigzag edge reconstruction modifies, not only the 
hoppings, but also the areas of the pentagons, heptagons 



In Fig. 11 ^a), we display the spectrum of a zigzag 



ribbon with real zz{57) edges in the presence of a per- 
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pendicular magnetic field. In contrast with the previous 
case, the two zero-energy levels are now split in energy 
and dispersive, crossing each other at the F-pointP^ 




Figure 11: Panel (a) shows the tight-binding energy spectrum 
of a zigzag graphene nanoribbon with real zz{57) edges with 
a width of 214 A (or N = 100 zigzag rows) in the presence of 
a perpendicular magnetic field, B = SOT. Panels (b)-(p) show 
the squared modulus of the lowest-energy levels wave func- 
tions for different values of ka [blue and red levels in panel 
(a)]. The green dashed vertical lines in (a), indicate the dif- 
ferent values of ka for which the edge states wave functions 
squared modulus were plotted in {b)-{p). 

The plots of the wave functions suggest a clear inter- 
pretation of this result. In graphene there is a bulk zero 
energy LL which cannot be affected by BCs, because the 
corresponding wave functions are localized in the bulk 
and do not reach the edges. And, in fact, one can see in 
Fig. pTj^g) regions of k with a flat energy level at zero 
energy; the plots of the corresponding wave functions 
(-0.9 < ka < -0.5 and 1.2 < ka < 1.5) show one lo- 
calized state inside the ribbon. The remaining states are 
edge states localized at its boundaries. In a real recon- 
structed edge, these states are dispersive in zero magnetic 
field, as we have seen in the previous section, and remain 
dispersive in a magnetic field: hence the lifting of the 
degeneracy and the level crossing at the F point, which 
involves states localized at opposite edges. On the other 
hand, in the simplified zz{57) ribbon, the edge states oc- 



cur at zero energy, as we have also seen. So the doubly 
degenerate zero energy state is either a zero energy bulk 
LL and an edge state, or two edge states, located at op- 
posite ends of the ribbon. This is confirmed by the plots 
of the wave functions. 

We now proceed to indicate briefly how these results 
arise from the Peierls substitution. We begin by consid- 
ering the appearance of a zero energy bulk Landau level 
(BLL). The recurrence relations for the amplitudes now 
involve matrices that depend on the row index, 

A(A:;n + l) = WA{n)A{k',n), (29a) 

B{k',n^l) = W^\n^l)B{k;n). (29b) 

Recall that A{k;n) and B{k;n) are notations for the 

column vectors, A{k;n) = [Ai{k;n),A2{k;n)] and 

B{k]n) = [Bi{k]n), B2{k]n)] ; the matrices I^a(^) 

and Wsin) are written in Appendix |c[ As before, these 
are commuting matrices, and have the common basis 
{u+,u_}; the eigenvalues, however, depend on the row 
index. 



^+(r) = -2e^^^/2 



cos 



ka 

Y 



^^{r) = 2ie^^^/2g- 



sm 



ka 

Y 



(r + l)7r- 



(r + l)7r- 



(30a) 

= (g(o)*. 

(30b) 



We can then rewrite Eqs. (29), for r?i > 2, as 



A{k;n) = E\{n^m)a-^{k;m)u~^ 

+ E^{n^m)a-{k;m)u~ ^ 

B{k;n) = E'^{n^m)l3-^{k;m)u~^ 

-f E^{n^m)f3-{k;m)u~ , 



(31a) 



(31b) 



where n > m^ a^ and (3^ are undetermined coefficients, 
while the quantities E^{n^m) and E^{n^m) are a short- 
hand for 



n-l 



^Ain,ni) = UeAir), 



E%{n,m) 



n 



Vi^bW 



(32a) 



(32b) 



As a function of the row index n, S^/^Jn,™) goes 



through a maximum when 



^Air) ^%{r) 



decreases 



below 1. These maxima are repeated periodically when 
n changes by 2n0, where n^ = (/)o/06 is the number of 
hexagons required for a total flux equal to a flux quan- 
tum. These multiple maxima are related to commen- 
surability effects between the lattice parameter and the 
cyclotron radius and are only important for unrealisti- 
cally high fields.^ For achievable values of the magnetic 
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field, n^ is mucii larger than the ribbon width, A/", (for 
B = SOT, 2n0 ~ 2000), and at most one maximum of 
'^A(B)i'^^ TTi) is located inside the ribbon, as shown in Fig. 
Assume, for instance, that that is the case for S^ at 



12 



riB- 



nA+ 



riA- 



ka 

2^ 

ka 

2^^ 
ka 
2^ 
ka 



ricj, 



Ucj) 



+ ^ h^0, 



q ]n^ = nB 



nB+ = —n^ 



6 ' 

2n0 
3 ' 

- +g I n^ = nB- + Y' 



1 <CnB- <A^- 1, 

(33a) 

(33b) 



q n^ = riB- 



(33c) 
(33d) 



where g is an integer. From Eqs. (33), we conclude that 



for reasonable values of the magnetic field and ribbon 
widths, at most, only one of the components will have a 
maximum inside the ribbon (of width N = 100). See, as 
an example. Fig. 12 
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Figure 12: Plot of |S^(n)| and |S^(n)| (that can be inter- 
preted as the amplitudes of the four components of the wave 
function in the proper basis of the matrices VK), given by 
Eqs. (32). The above quantities were plotted for B = SOT 
and ka — —0.716. In the (a) panel, the \Ej^.^(n)\ are nor- 
malized over the region n G [—500, 500], while in the {h) panel 
they are normalized over n G [0, 100]. The red box in panel 
(a) signals the region where n G [0, 100]. 

Moreover, the amplitude P-{k^n) will decay exponen- 
tially to very small values at the edges; to exponential 
accuracy, the BCs, whatever they may be, are trivially 
satisfied by choosing a^ = «_ = /3+ = 0; this then is a 
BLL, where the wave function exists only in one of the 
sub-lattices and is localized away from the edges. These 
BLLs occur irrespective of the type of edge. However, 
when k changes and the LL center approaches the edge, 
the BCs come into play, differentiating the various situ- 
ations. 



Let us now consider the appearance of the edge states 
in these results. The general BC for a reconstructed 
zigzag edge with SW defects at the n = end, may be 
written as a(/c;2) = A1/3(/c;2), where Ai is defined in 
Appendix [Cj Eq. ( |C4[ ), with an analogous expression for 
the edge in n = TV - 1, /3(fc; TV - 3) = M'a{k; N - 3). 

We will start by assuming that the ribbon is terminated 
with a simplified zz{57) edge (hi = v = h[ = 1 and (j)j = 
(^5 = 2(j)(iQ = 206)- In such a situation, we have A^++ = 

M. =0, a result which uncouples the components «+ 

and j3- from a- and /3+ 

a+(A:;2) = A^+_/3_(A:; 2), (34a) 

a_(/c;2) = A^_+/3+(/c;2). (34b) 

As a consequence, every time we have a zero energy 
BLL (living away from the edges), we will also have one 
other solution of zero energy, now localized at the edge. 
Let us take as an example, the case where ka = —0.716, 
for which the S^,^ are depicted in Fig. (12). To expo- 
nential accuracy, the BCs involving /3_ anaa+ are triv- 
ially satisfied at both edges choosing a+ = 0. Those in- 
volving /3+ and a- , are satisfied at the upper edge choos- 
ing a_(/c; 2) = M. h/^+(^5 2), being satisfied at the lower 

edge to exponential accuracy. In the real space [see Eqs. 
(31)], we will have a BLL localized only on the B sub- 
lattice and an edge state around the edge at n = 0, living 
in both sub-lattices with different localization lengths. 

When the value of ka is increased, the maxima of S^ ,^ 
move to higher values of n. At a certain point, the max- 
imum of S^ will be such that ub- > N — 1, and then 
there will be no maxima inside the ribbon. In such a case, 
the maxima of S^ and S^ closer to the ribbon, will be at 
n > N —1^ while the maxima of S^ and S^ closer to the 
ribbon, will be at n < 0. In this case, the BCs involving 
P- and a+ will be satisfied at the lower edge choosing 
l3-{k; N -3) = M'_^a^{k; TV - 3). At the upper edge, 
the BC will be obeyed to exponential accuracy. The con- 
verse needs to be done regarding the BCs involving /3+ 
and a-. Consequently, for —0.5 < ka < 1.2, there will 
be zero-energy solutions localized at both edges, living 
in both sub-lattices with distinct localization lengths in 
each sub-lattice. 

If, on the contrary, we start decreasing the value of ka 
from ka ~ —0.72, the maxima of S^^^ moves to lower 

values of n, and at a certain point, the maximum of S^ 
will be such that ub- < 0. In such case, the maxima of 
S^ and S^ closer to the ribbon, will be at n < 0, while 
the maxima of S^ and E\ closer to the ribbon, will be 
at n > TV — 1. In this case, it will not be possible to 
satisfy the BCs non-trivially and consequently, there will 
be no zero-energy solutions in this region, as can be seen 
in Fig. 
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When instead of simplified zz{57) edges, the ribbon 
is terminated with real zz{57) edges, the matrix A^ is 

modified, and A^++ 7^ A^ 7^ O7 resulting in a BC 

coupling all the components a± and P± 
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a+(A:;2) = M++/3+(A:; 2) + A^+_/3_(A:; 2), (35a) 

a_(fc;2) = A^_+/3+(fc;2)+;W__/3_(/c;2). (35b) 

To grasp the implications of this modification, con- 
sider for instance the case where ka = —0.716, depicted 
in Fig. [T2J where a BLL is present in the f3- mode; since 
/3_(2) ^ 0, the BC imply all three remaining amplitudes, 
Qf+, a_and /3+ to be non zero, if there is to be an edge 
state in addition to the BLL. But the a+ mode grows 
as n increases, whereas the other two decrease; as a re- 
sult the BCs will be violated at the opposing edge. In 
conclusion, BCs can no longer be satisfied with zero en- 
ergy edge states, which become dispersive, whereas zero 
energy BLL still occur. This accounts for the lack of 
zero energy doubly degenerate state in ribbons with real 
reconstructed zz{57) edges. 



III. CONCLUSION 

We have discussed in detail the effect of edge recon- 
struction on the characteristics of low energy edge states 
in graphene ribbons. In the case of Stone- Wales zz{57) 
reconstructed zigzag edges, we find a new type of edge 
state originating from the doubling of the unit cell along 
the edge, brought about by the edge reconstruction. This 
new type of edge state has the following features: (i) the 
states are in general dispersive, although specific values 
of the tight-binding model parameters allow zero energy 
states; (ii) the wave- function, even for the semi-infinite 
ribbon, has non-zero amplitudes on both sub-lattices; 
(iii) close to the Dirac points, the wave function ampli- 
tudes have two components decreasing with the distance 
from edge with different decay lengths, one of which re- 
mains finite, of the order of the lattice parameter, even 
at the Dirac point, while the other diverges. The disper- 
sion of the edge states should lead to a charge transfer 
between bulk and edges (self-doping), which, for realistic 
values of the tight-binding parameters, leaves the edges 
negatively charged. 

In the presence of a magnetic field, one still finds zero 
energy bulk Landau Levels, as was to be expected, since 
these are insensitive to the edges; however, in contrast 
to pristine zigzag ribbons, where the zero energy LL is 
degenerate with an edge state, this in no longer true in 
ribbons with reconstructed edges, since the edge states 
are, in general, dispersive. 
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Appendix A: Tight Binding equations and boundary 
conditions 



In this appendix, we write the tight-binding equations 
for the amplitudes at the sites near one edge, n = 0; these 
will determine the boundary conditions (BCs) that must 
be satisfied by the bulk solutions. For clarity, we begin 
by considering zero energy states. We will argue that the 
BCs adequate for low energy states, \e/t\ <C 1, are the 
same as for zero energy states. 

The tight-binding equations at the sites of Ai{m^O) 
and A2{m^0) have the form: 

h2A2{m-0)^hiBi{m;0) = 0, (Ala) 

h2Ai{m,0)^hiB2{m;0) = 0. (Alb) 

It will be useful to express these in matrix form; after 
Fourier transforming in the m index, {k is the wave vector 
along the edge). 



A{k;0) 



ho 



a,B(fe;0), 



(A2) 



where o-^ is a Pauli matrix. For the Bi{m;0)^ B2{m]0) 
sites, 

h^B2(m - 1; 0) + /iiAi(m; 0) + vAi(m] 1) = 0, (A3a) 
h^Bi{m + 1; 0) + hiA2(m] 0) + vA2(m] 1) = 0. (A3b) 

Using Bloch's theorem, we can cast this in the form 

hi 



(A4) 



A(fc;l) + - 


-^A(fc;0) 

V 






/i4 

V 


g-2i/ca Q 

e^'^"" 


(T^B(fc;0)=0. 


Using Eq. (|A2 ) in this one, 


A(/c;l)+7^c^^B(A;;0) = 0, 


where 




" h' 


l-h2h4e-^'''^ 




h2V 




7^:=- 




h2V 

^ 


: 



(A5) 



(A6) 



is a matrix that depends on k. 

With a similar procedure for the sites Ai{m;l), 
A2{m; 1), Bi{m; 1) and B2{m^ 1), we obtain 



A(fc;2) 
B(fc;0) 



W^A(fc; 1), 
WBB(fc;l), 



with 



Wa 



Wb 



h[ h[ 


■) 


' h[ h'^e-^^'"'' 


h[ h'. 





(A7a) 
(A7b) 



(A8a) 
(A8b) 
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using Eqs. (A5), we arrive at 



Naturally, this reduces to Eq. (A13a) if the right hand 



A{k; 2) + >VA7^c^^>VBB(^; 1) = 0. (A9) 

Beyond the first row (n > 1), it is simple to get 

A(/c;n + l) = WAA{k,n), (AlOa) 

B(A:;n + l) = W^^B{k,n), (AlOb) 

where 

Wa = - J,^ \ , (Alia) 

1 p — 2ika 

Wb = - \\ , (Allb) 

In summary, after Fourier transforming in the m vari- 
able, the tight-binding equations for a semi-infinite rib- 
bon with zz{57) reconstruction are (n > 1) 

A(fc;2) = -WA7^a^WBI^BB(fe;2), (A12a) 
A(/c;n + l) = WAA{k,n), (A12b) 

B(A:;n + l) = W^^B{k,n), (A12c) 

The last two are the bulk recursion relations, while the 
first one contains the BC. 

We now generalize these equations for states of finite, 
but low, energy. We argue that only the bulk equations 
are changed, the BCs remain the same, i.e., 

A(/c; 2) = -WA7^a^Wi?T^BB(fc; 2), (A13a) 

A(fc, n + 1) - WA^ik, n) = - (-\ B(/c, n + 1), (A13b) 

B(fc, n) - WsBik, n + 1) = - (-) A(k, n). (A13c) 

Let us put back the energy in the equations for the am- 
plitudes near the edge, 

/i2A2(m;0) + /ii5i(m;0) = - (^) Ai(m;0), (A14a) 
/i2Ai(m,0) + /ii52(m;0) = - (^) A2(m;0), (A14b) 



so Eq.(A2) becomes. 



A(m; 0) + T;^a^B(m; 0) = — (--\ a^A(m; 0). (A15) 

112 tl2 \ t / 

This shows the pattern that we have to repeat in Eqs. 
(TASI) through to Eqs. (TATI). Instead of Eq. (|A13a|), we 



obtain, 



A(fe; 2)+WA7^a,WBl^BB(fc; 2) 

-i-d 



-WAB{k;0) - ^WACTa:A{k;0) 
V h2V 



-WA7^a^A(A:;l) + B(A:,l) 



side is set to zero. The important point is that, for the 
values of the parameters listed in Table |Tj the matrix 
Wa'^ct^Wb^b has one finite eigenvalue in the entire 
range of /c, whose modulus is always larger than about 
1.3. This means that, to lowest order in (— e/t), we are 
justified in neglecting the RHS of this equation, and use 
the same BC as for zero energy states. This is a valid 
approximation for states with \e/t\ <C 1. 

Now we change basis to rewrite these equations in the 
eigenbasis of Wa and Wb^ [see Eqs. (14)] 



1 


^-ika 




V2 


1 


7 


1 


_^—ika 


V2 


1 





(A17a) 
(A17b) 



The coordinate transformation is defined by the matrix 
U given by 



U- 



V2 



_^ika I 



(A18) 



The BC in the new basis, becomes 



(x{k;2) = -UWAna^WBWBU'^^{k;2) 

= M(k)f3{k-2). (A19) 

and the bulk equations, 

a^{k] n + 1) - iA^cj{k] n) = - (-] (3a{k; n + 1), 

(A20a) 

PAk; n) - {^iyp{k- n + 1) = - Q) a^k- n). (A20b) 

The matrix M{k) can be calculated explicitly, since 
all the matrices intervening in its definition were given 
above, but its long expression is not particularly enlight- 
ening. 

Appendix B: The loAv-energy edge states 

We now sketch the calculation of th e low energy edge 
states for the problem set by Eqs.(A13) in a semi-infinite 



(A16) 



ribbon. For solutions that decay away from the edge, 

a,{k; n) = e^«-("-2)a^(fc; q,), (Bla) 

/3,(fc;n)=e^«-("-2)/3,(fc;q,), (Bib) 

the equations for the amplitudes in the bulk become 

(1 _ e-'-i-a)a,{k; q^ = - (^) Mk; Qa), (B2a) 
(1 - e^''-(a)*)/?.(fc; in = - Q Mk; q.)- (B2b) 
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The energy must be given by 

2 



(!) =(l-e-^«-C5)(l-e+*''(ar) 



(B3) 



Expanding the RHS, and given the fact that the energy 
must be real, we conclude that ^ [e~*^^'^^^] = 0, which 
is equivalent to e"**^^^^ = ± |^^| e^*^"^. This allows us to 
rewrite Eq. (B3) as 



(^)' = l + leSl'T2|eS|cosh(9q-). (B4) 

Low energy solutions, with \e/t\ <C 1, correspond to the 
choice of the minus sign in this expression. From this, we 
can write the energy expression as 



-(l-iaie^''") 



1 



(B5) 



On the other hand, the energy can be eliminated from 
Eqs. (B2) to obtain. 



3^g. 



''n\ 



1 _ e-a?^ 1^2 



aa{k,qcr) 



(B6) 



This result shows that the values of ^^cr are determined 
if we fix the amplitude ratios, s^r, i-e., if we take as BCs 
for the two a = +, — , chains 

To determine the value of the energy we use the fol- 
lowing conditions: (i) the values of 5+ and 5_ are related 
by the BCs [Eq. (|A1^], 



M 



++ 



■det[A^]5_ 



1-M- 



(B7) 



(ii) their values must be such that the RHS of Eq. ( B6 ) 
is independent of a. Hence, we determine ^(7_and ^(7+, 
as a function of 5_ (using the value ofs+ given by Eq. 
(B7), calculate the energies from Eq. (B5) for a = +, — , 
and vary 5_ until the two energies match; as long as 
jg^g'^ I < 1, this constitutes the solution of our problem. 
Note that the sign of the ener gy, is determined by the 
hopping amplitudes trough Eq. (B7). The BCs we used 



are only valid for \e/t\ <C 1. As a consequence, we can 
expect that this analytical construction of edge states will 
only be valid near the Dirac points {ka = ±7r/3), where 
this condition is fulfilled. 



Appendix C: Recurrence matrices Avith magnetic 
field 



When a perpendicular magnetic field is applied perpen- 
dicularly to the ribbon, in the bulk, the matrices Wa(^) 



and WsiTi) read 



WA{n) = 



«("+i)-|§ 



i{n+l)nl 



i(2fea-(n+l)^|fi) i{n+l)ni 



(Cla) 



Wein) 



*(»+l)-|t g-<2fea-(„+l),r|a) 









(Clb) 



where ^6 is the magnetic fiux through an undistorted 
hexagon. 

Moreover, the matrices around the upper edge, W^, 
W^ and 7^, are given by 



w^ 



m 



hi e '^o h\ e 



4>o 



7 / i7T-^ 



n = - 



vh2 







vh2 



00 



(C2a) 

,(C2b) 

(C2c) 



where = 2ka — tt^ — tt^ and ^7^, (j)^ and ^^g are the 
fiuxes of the magnetic field across the upper heptagons, 
pentagons and distorted hexagons (see Fig. |9|, while (j)o 
is the fiux quantum. The matrix associated with the 
boundary at the upper edge, o^, reads 



CTx = 



u 



. 0- 
e '^o 

^0 



(C3) 



If we take the energy to be zero, and change to the 
proper basis, the BC for the edge at n = becomes 

a(fe;2) = -UW'^n^W^WB{2)U^^{k;2) 

= M{k)^{k;2). (C4) 

The proper basis of matrices Wa and Wb is {u~^ ^ u~} 
defined in AppendixjAj In the proper basis, the equations 
for the bulk amplitudes, read 

- a^(/c;n + 1) +^^(n)a^(/c;n) = - f - j /3^(/c; n), 

(C5a) 
-^{n)f3^{k;n)^f3^{k;n-l) = -Q)a^(A:;n), 

(C5b) 



where the ^^/^ are defined in Eqs. (30) 
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The numerical diagonalization of the tight-binding Hamil- 
tonian was performed using the tools of LAPACK numer- 
ical library. 

The alternative possibility for zero energy states in the 
range where |^^| > 1, and a+(/c) = a-(k) = 0, requires 
det [A^(A;)] = 0; we found no relevant limits where this is 
verified. 

We have also confirmed numerically the prediction that 
edges states have zero energy when /i? — /i2/i4 = 0, though 
we do not present the corresponding data. 
The whole spectrum is shifted in ka, because we have set 
n = at the upper edge (where n is the label of the zigzag 
rows). If we have set n = to the center of the ribbon, the 
shift would disappear.'^ 



